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1.  Introduction 


It  is  well  known  that  multipath  fading  can  significantly  effect  the  link  reliability  in  a 
communications  system  or  target  detectability  in  the  case  of  radar.  The  path  between 
a  transmitted  and  a  receiver  is  often  obstructed  by  natrual  or  man-made  obstacles 
such  as  hills,  buildings,  atmospheric  layers,  trees,  rain,  fog,  etc.  Propagation  outage 
due  to  multipath  fading  depends  in  a  complicated  manner  on  propagation  climate, 
terrain  features,  path  length,  radio  frequency,  and  fade  margin.  In  the  case  of  atmo¬ 
spheric  multipath  fading,  interference  due  to  two  or  more  super- refracted  rays  arriving 
at  the  receiver  via  different  paths  can  lead  to  a  complete  loss  of  signal.  Other  patho¬ 
logical  phenomenon  such  as  obstruction  fading  (caused  by  sub-refractive  atmospheric 
effects)  and  ducting  (caused  by  extreme  super-refractive  effects)  are  also  possible. 
Such  phenomenon  are  more  common  in  wrarm  and  tropical  climates,  particularly  near 
shores,  w'here  elevated  inversions  are  formed  easily  due  to  the  large  temperature  and 
partial  pressure  differentials.  Reflection  multipath  fading,  which  is  due  to  interfer¬ 
ence  between  the  direct  and  the  ground  reflected  ray  depends  strongly  on  the  terrain 
geometry  and  ground  constants.  Moreover,  elevated  terrain  features  could  completely 
mask  a  receiver  from  a  transmitter  leading  to  severe  loss  of  signal  (in  some  cases,  it 
is  advantageous  to  site  antennas  behind  hills  to  provide  shielding  against  undesirable 
interference).  It  is  very  important  to  assess  the  effects  of  environment  on  the  link.  A 
computer  model  that  can  take  into  account  a  given  refractive  index  profile,  terrain 
elevation  data,  and  varying  ground  parameters  will  be  very  helpful  in  predicting  the 
link  performance. 

In  this  report  we  present  formulation  details  for  an  efficient  numerical  solution  of 
wave  propagation  in  an  inhomogeneous  atmosphere  and  over  irregular  terrain  using 
parabolic  equation. 

Unlike  all  other  previous  formulations  of  the  parabolic  equation,  we  will  use  a  modified 
Helmholz  equation  for  propagation  in  an  inhomogeneous  atmosphere  as  suggested  by 
Maxwell’s  equations  (all  previous  formulations  use  a  Helmholtz  equation  which  is 
only  true  for  fields  in  a  homogeneous  medium).  Because  the  parabolic  equation  is  a 
full- wave  method,  it  will  include  all  aspects  of  wave  propagation  such  as  reflection, 
refraction,  diffraction,  and  surface  wave  propagation.  In  this  respect  it  is  far  superior 
to  the  commonly  used  ray  method. 

Parabolic  equation  approximation  to  an  elliptic  partial  differential  equation,  which 
the  true  fields  satisfy,  has  proven  to  be  a  viable  approach  for  studying  propagation 
problems  in  underwater  acoustics.  The  method  is  just  gaining  popularity  with  the 
electromagnetic  community.  Although  the  parabolic  equation  regards  waves  as  es¬ 
sentially  traveling  one-way,  it  allows  a  rapid  solution  of  the  fields  by  way  of  marching 
along  the  range  starting  from  an  initial  range.  Another  advantage  of  the  PE  method 
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compared  to  the  ray  methods  is  that  it  is  valid  even  in  the  shadow  region  where  the 
simple  ray  methods  completely  break  down.  Furthermore,  it  appears  to  be  the  only 
practical  method  for  predicting  propagation  over  long  ranges  (greater  than  1  km)  over 
a  wideband  from  HF  (a  few  MHz)  through  SHF  (a  few  tens  of  GHz).  The  method  is, 
however,  not  without  limitations.  In  its  standard  form,  the  accuracy  of  the  method 
is  limited  to  waves  traveling  essentially  within  ilO  from  horizontal,  furthermore, 
treatment  of  the  boundary  conditions  on  the  uneven  terrain  is  difficult. 

The  method  we  propose  to  use  will  attempt  to  remove  both  of  these  deficiencies. 
Firstly,  we  use  a  Helmholtz-like  elliptic  equation  to  describe  the  fields  and  arrive 
at  a  wide-angle  parabolic  equation  subject  to  certain  approximations.  To  facilitate 
imposition  of  the  boundary  conditions  on  the  irregular  terrain,  the  equations  will  be 
transformed  to  a  body-fitted  curvilinear  coordinate  system.  The  PE  will  be  solved 
using  finite  differences  on  a  non-rectangular  mesh.  This  is  a  major  departure  from 
previous  approaches  which  will  not  only  make  the  method  more  efficient  but  also 
more  accurate. 

In  Section  2,  we  derive  the  exact  equations  satisfied  by  2D  fields  in  an  inhomogeneous 
atmosphere.  Impedance  boundary  conditions  are  used  to  characterize  the  ground.  In 
Section  3,  we  present  details  on  the  impedance  boundary  conditions.  Starting  from 
the  exact  equations  presented  in  Section  2  for  the  fields,  we  derive,  in  Section  4,  a 
parabolic  equation  (PE)  valid  for  narrow  angle  propagation.  This  case  is  termed  as 
the  standard  PE.  The  standard  PE  is  generally  valid  for  propagation  angles  that 
are  within  ±10°  from  horizontal.  To  accomodate  waves  traveling  at  larger  angles,  we 
present  the  derivation  of  a  wide  angle  PE  in  Section  5.  To  truncate  the  computational 
domain  we  derive  boundary  conditions  on  an  upper  boundary,  which  are  termed  as 
the  tropospheric  boundary  conditions.  The  derivation  is  based  on  the  solution  of 
Schrodinger  type  parabolic  equation  in  a  quarter  plane  x  >  0,  y  >  0.  This  is  presented 
in  Section  6. 

For  an  efficient  numerical  implementation,  we  transform  the  differential  equation  and 
boundary  conditions  to  a  curvilinear  coordinate  system.  This  is  presented  in  Section 
7.  Details  on  the  numerical  generation  of  the  curvilinear  coordinate  system  are  given 
in  Section  8.  Finally,  in  Section  9,  we  present  steps  leading  to  a  finite  difference 
solution  of  the  equations  using  a  Crank-Nicholson  implicit  scheme. 
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2.  Solution  of  2D  Fields  in  an  Inhomogeneous  Medium 


Consider  an  electric  source  producing  fields  in  an  inhomogeneous  region  as  shown 
in  the  figure  below.  Let  us  assume  that  both  the  sources  and  the  medium  are  two- 
dimensional  in  nature  in  that  all  quantities  are  independent  of  the  2-coordinate.  As 
in  the  case  of  a  homogeneous  medium  the  fields  can  be  decomposed  into  a  TEZ  case 
(vertical  polarization)  and  a  TMZ  case  (horizontal  polarization).  It  is  assumed  that 
propagation  takes  place  in  the  xy-plane. 


From  Maxwell’s  equations,  we  have 

V  x  E  =  -jufiH  (1) 

V  x  H  =  jutE  +  J  (2) 

In  the  case  of  vertical  polarization,  the  fields  could  be  written  in  terms  of  the  z- 
component  of  the  magnetic  field,  H  =  zHz ,  (TEZ  fields).  Substituting  into  (2)  we 
have 


V  x  H  =  V  x  (. zHz )  =  V//2  x  z  =  jutE  +  J  =>  eE  = 


V//2  x  z  -  J 


In  a  source-free  environment,  we  have 


eE  = 


VHZ  x  z 


Substituting  into  (1)  we  get 


1  „  /V//2  z  _  fVH.\  ft 

Vx£  =  —V  x  - x  z )  =  -—V  •  - )  =  -jujfiH 

ju)  \  e  J  ju  \  e  / 
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The  equation  satisfied  by  the  magnetic  field  is  then 


V 


LJ  = 


+  u>2[iHz  =  0  Vertical  Polarization 


Once  the  magnetic  field  is  determined  the  electric  field  is  given  by 

"I  z  x  VHZ 


Note  that  Hz  does  not  satisfy  the  Helmholtz  equation  unless  e  is  constant. 
For  horizontal  polarization  on  a  similar  analysis  with  E  =  zEz  shows  that 


y  /  je  \ 

- i  -)-  u2eEz  —  0  Horizontal  Polarization 


If  the  medium  is  non-magnetic,  fi  =  p0  and  Ez  satisfies  the  Helmholtz  equation.  We 
will  characterize  the  ground  in  terms  of  impedance  boundary  conditions  [3]. 

We  may  combine  the  vertical  and  horizontal  cases  shown  in  (3)  and  (5)  into  an 
equation  of  the  form 

V  •  (aVt/O  +  /?</>  =  0  (6) 


where 


TE  or  Vertical  Pol. 


TM  or  Horizontal  Pol. 


j3  =  akln2(x,  y)  >  0 


*l>  = 


Hr  TE  Pol. 


Ez  TM  Pol. 


The  quantity  n(x,y )  is  a  position  dependent  refractive  index  of  the  medium.  The 
partial  differential  equation  (PDE)  given  in  (6)  is  elliptic,  for  we  have 

f  dip\  d  (  dip\ 


dx)  +  dy  \Q  dy 


dx 


5 


+  =  0 


or 

'dV  d2i!>  1  d x(>  da  di>  da  _ 

a[d^  +  w\  +  fo^^^^ 

The  discriminant  for  the  above  PDE  is  -1  <  0  implying  elliptic  nature  of  the  equation. 
Equation  (6)  could  be  expressed  as 

„  1  da  ,  1  da  (3  .  .  nN 

VV+-n-^+-^*+^  =  0  (l0) 

a  dx  a  ay  a 
For  a  non-magnetic,  loss-less  medium,  we  have 


- ,  Vertical  Pol. 

tQtr 


Horizontal  Pol. 


so  that  for  vertical  polarization 


i  da  _  £  d_n 

a  dx  r  dx  \er 

2  dn 
n  dx 


\  1  dtr  _  1  dn 2 

)  tr  dx  n 2  dx 

(■ n  is  the  refractive  index) 


Similarly, 


Letting 


ai(x,y) 


1  da  2  dn 

a  dy  n  dy 


2  dn 
n  dx 


Vertical  Pol. 


Horizontal  Pol. 


a2(x:y)  = 


2  dn 
n  dy 


- —  Vertical  Pol. 


Horizontal  Pol. 


equation  (10)  may  be  expressed  in  the  form 


V2^  +  ai%l)x  +  a2^y  +  kln2ij)  =  0 
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3.  Impedance  Boundary  Condition 

Impedance  boundary  condition  relates  the  tangential  components  of  electric  and 
magnetic  fields  at  the  interface  of  two  media.  If  £  is  a  unit  normal  and  s  is  a  unit 
tangent  as  shown  in  the  figure  below,  the  boundary  conditions  are  given  by  [3] 

s  =  x  cos  6  +  y  sin  9 
v  =  —  x  sin  9  +  y  cos  9 

x  —  s  cos  9  —  v  sin  9 
y  =  s  cos  6  +  u  cos  6 

(12) 

where  As  =  Zs/rj0  is  the  surface  impedance  normalized  to  the  free  space  value  rj0. 
The  equation  may  also  be  expressed  as 

(0  ■  E)v  -  E  =  -t)0As{)  x  H  =>  v  x  E  =  770 Asi>  x  (i>  x  H) 

=  VoAs  [(u  •  H)v  -  H\  (13) 

The  surface  impedance  is  determined  from  the  intrinsic  impedance  of  the  medium 
by  considering  plane  wave  reflections  from  the  interface.  The  complex  propagation 
constants,  71,  72,  and  the  intrinsic  impedances  r)i  and  rj2  in  terms  of  the  media 
constants  are  indicated  in  the  figure. 


7i  =  juHrif*o(<7i  +  jue0eTi) 


V  2  -  Vo 

According  to  Snell’s  law,  we  have  71  sin  6,  =  72  sin  6t 
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The  plane  wave  reflection  coefficients  for  the  vertical  and  horizontal  polarizations,  Ry 
and  Rh  are  [7] 

77o  SeC  Of  TJi  SeC  0 ,  n  r  T  J  *7//  f\ 

Rh  =  -A - l - Li - -  =>  Surface  Impedance  Zs  -  r)2  sec  0t 

r/2  sec  0t  +  7/1  sec  0,- 


T]2  COS  —  TJ\  cos 

T\i  COS  —  7/1  COS  Oi 


H  T)2  secOt  _  7?2 


Af  - 


ZVS  =  7/2  COS  0t 


Vi  V\ 


=  !_  —  sin  # 


j  Pt2^tc1 
P rl  ^rc2 


.  PrRrcX  2  / 
1 - COS  xp , 

Pr2^rc2 


.1/  I  Pt2^tc1  i  PtX^tcI  2  / 
A^  =  - 1 - “COS  xp  i 

V  Prlerc2  Pr2Ec2 


For  the  special  case  of  p-y  =  p2  =  ^o,  —  0; 


A?  = 


A.v  = 


■  1- 


er2  “  Jcrr2 


•  cos2  xp. 


Cr-2  ~  j&T2 


Cr2  “  7^2 

i  2  / 

1 - ; - COS  Xpi 

Ct-2  “  7^2 


For  normal  incidence  xpi  =  90° 


=  Av 

S  s 


tr2  ~  J(Jr2 


For  the  2-D  case,  impedance  boundary  conditions  for  vertical  polarization  can  be 
simplified  as 

v  x  E  —  7/oA^  (z>  •  H)v  -  #]  —  -rjoA^H 

Taking  a  dot  product  on  both  sides  with  z 
z  -{v  x  E)  =  —t)0A^Hz 

Since  2  x  r  =  — s,  we  have 
-s  ■  E  -  -rj0A]'Hz 

— * 

Substituting  from  equation  (4)  for  E ,  we  get 

,  zxVHz 

s  - - : - =  -rj0AsHz 

jue 
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v  •  V#2  =  juerjoAl  Hz  =  jk0erA^Hz 
i.e.,  the  above  can  be  put  conveniently  in  the  form 

IBC  Vertical  Pol. 

A  similar  analysis  for  horizontal  polarization  yields 

IBC  Horizontal  Pol. 

This  may  also  be  obtained  by  resorting  to  duality. 

For  a  perfectly  conducting  material,  we  have  =  0  and 

—  0  Vertical  Pol. 

dv 

Ez  —  0  Horizontal  Pol. 


(16) 

(17) 
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4.  Standard  PE  Derivation 


We  will  make  some  approximations  and  cast  (11)  in  the  form  of  a  parabolic  equa¬ 
tion  which  permits  a  rapid  numerical  solution. 


Let  i>(x,y)  =  — j=-u(x,y) 
v  x 


Then 


4>x  =  ( ux-jk0u 


u  \  e 


2x3/2  J  yfi 


{ux  -  jk0u )- 


x  — >  oo 


Tpy  =  ul 


\fx 


—  U xy  ~  (U-xy  jkfQUy^' 


vyy  ~  ayy  ^  ’ 
Substituting  into  (11)  we  get 


(uxx  2  jkoUx  k0u)- 


uxx  +  Uyy  -  2jk0ux  +  a\Ux  -  2jk0axu  +  a2uy  +  (n2  -  1)/cqU  -  0 

°r  /  fll\ 

Uxx  +  Uyy  +  (fli  -  2jk0)ux  +  a2uy  +  [n2  -  1  -  2j— J  Ku  =  0 

If  now  we  impose  the  approximation  that 


we  obtain 


<  (a2  +  4 k2)1/2  |ur|  , 


(ax  -  2jk0) 


Uyy  +  a2Uy  T  2  -  1  2 j  ^  U 


u(x)  = 


(2/cq  +  ja  i ) 


^  d 2 

Gl  +  a2^  +  V 


This  is  the  exact  form  of  narrow  angle  PE  approximation.  We  would  also  like  to 
express  the  impedance  boundary  condition  in  terms  of  the  u  functions. 


x  =  s  cos  0  —  v  sin  6 
dx 

xv  =  ^~  =  v-Vx  =  v-x 

ov 

—  —  sin  6 

v  =  —  £  sin  0 -f  $  cos  0 
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For  vertical  polarization  we  have 

dHz 


du 


jk0tTAs  Hz  =  0 


dHz 

du 


jk0n2AvsHz  =  0 


-jkox 


Hz(x,y)  =  ~=-u{x,y) 
=  (u„—j i  kQxvu 


ux,,  \  e 


,-jk ox 


du 


(u„  -  jkQxuu) 


2x3/2y 


-jkox 


x  — ►  oo 


-  jk0(n2A ]  +  x„)u  =  0 


IBC  Vert.  Pol. 


Similarly  for  horizontal  polarization,  we  have 


u„  -  jko  \^H+Xt'j  u  =  0 


IBC  Horz.  Pol. 


We  combine  the  two  by  defining 

-jk0{n2A^  -  sin0)  Vert. 


ci  =  1 


-jk0 


1 

A? 


sin  6  Horz. 


and  writing  as 


u„  +  Ci«  =  0 


(19) 


The  parabolic  equation  given  in  (18)  is  valid  for  propagating  angles  close  to  horizon¬ 
tal  (±10°  in  practice)  [1].  To  accomodate  waves  at  higher  angles  we  would  need  a 
wide  angle  parabolic  equation  whose  derivation  is  accomplished  through  a  pseudo¬ 
differential  operator  formalism  [1]. 
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5.  Wide  Angle  PE  Derivation 


Let  us  assume  that  the  media  constants  are  independent  of  horizontal  range  so  that 
a{x,y)  =  a(y),  n{x,y)  =  n(y).  We  then  have  from  (17) 

uxx  +  Uyy  -  j2k0ux  +  —  —  +  ( n2(y )  _  l)^ou  =  0 

Let 


Using  this  notation,  the  above  PDE  may  be  expressed  as 

P 2  -  2 jk0P  +  {Q2  -  l)fco]  u  =  0 

which  we  may  factorize  as 

(. P  -  jh  -  jk0Q)(P  -  jk0  +  jk0Q)u  =  0  (20) 

To  see  this  we  expand  the  operator  on  the  left  hand  side  to  get 

P 2  -  jk0P  -  jk0QP  -  jk0P  -k\-  klQ 
+  jk0PQ  +  klQ 2  +  k20Q 

Since  PQ  =  QP ,  we  have  the  desired  result.  The  first  operator  in  (20)  denotes 
an  incoming  wave  (w.r.t.  x)  and  the  second  an  outgoing  wave.  We  retain  only  the 
outgoing  wave  to  obtain 

Pu  =  -jk0{Q-l)u  (21) 

The  square  root  operator  Q  is  global  in  nature  and  we  would  like  to  make  some 
approximations  to  derive  a  local  operator  from  it.  (It  is  global  because  when  expanded 
in  terms  of  series,  it  will  contain  terms  of  all  derivatives).  Now 

r  2  II  da  d_  l_  (P] 1/2 
Q-[n+kladydy  +  kl  dy\ 


Id 2  1  da  d  2 

^  ~  \  k20  dy2  +  kl&  dy  dy+U 


Pseudo  differential 
operator  in  y 
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which  may  be  expressed  as 


Q  =  n/TTv7 


where 

2  1  da  d  1 

V  ~  H  1  +  k20a  dy  dy  ^  k20  dy 2 

=  a  small  operator! 


Treating  V  as  an  algebraic  factor  <  1,  we  may  derive  the  following  rational  approxi¬ 
mation  (pade(  1, 1)) 

Q  =  y/l  +  V  «  - r —  (Claerbout) 

1  +  4V 

4  +  3V  2V 

-  4  +  U  “  +  4  +  V 


so  that 

Q- 1 

Substituting  into  (21)  we  arrive  at 


2V 

4  + V 


Pu  =  j/c0(Q  -  l)u  PU  - 


—2jk0V 
4  +  V  1 


or 


Using  the  fact  that  a2 


(4  +  V)Pu  =  -2jk0Vu 

1  da 

-  — ,  we  get 

a  dy 


(n2  +  3)/cq  +  a2^-  + 


a2 

%2 


Ux  2  J  /uQ 


w2  5 

U*b  +  + 


dy 2 


u 


(22) 


The  above  equation  is  a  wide-angle  parabolic  equation  valid  for  propagation  up  to 
±20°  [1].  Other  approximations  could  be  obtained  by  considering  higher  order  pade 
approximants. 
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6.  Boundary  Condition  on  the  Upper  Boundary 

To  truncate  the  computational  domain,  we  consider  a  point  high  enough  where  the 
atmosphere  is  homogeneous  with  n  =  1.  The  governing  equation  in  the  homogeneous 
region  becomes 

j 

ux  —  2kQUyy 

Let  us  derive  boundary  conditions  on  a  horizontal  interface  y  =  y0.  For  the  sake 
of  simplicity  we  will  derive  boundary  conditions  of  the  mixed  type  on  the  interface 
y  =  0  (instead  of  y  =  y0)  given  the  initial  data  on  the  line  x  =  0  and  boundary  data 
on  the  line  y  -  0.  Our  derivation  is  based  on  the  use  of  Fourier  sine  transforms  as 
suggested  in  [2].  Although  the  basic  philosophy  of  our  approach  coincides  with  that 
in  [4],  some  of  the  details  and  the  final  results  are  slightly  different  from  the  latter. 

Consider  the  parabolic  equation  uyy  —  2 jkux  =  0  in  a  homogeneous  region  x  >  0, 
y  >  0,  where  A:  is  a  complex  constant, 


Ux 


x 

subject  to  the  initial  condition  u( 0,  y)  =  f{y),  0  <  y  <  oo,  and  the  boundary  condition 
u(x,0)  =  g{x),  0  <  x  <  oo.  We  assume  that  u(x,  oo)  — *  0  and  uy(x,oo)  — *  0.  The 
equation 

uyy-2jkux  (23) 

is  of  Schrodinger’s  type.  We  will  treat  the  lossless  case  having  a  real  value  of  k  as 
the  limiting  case  of  the  lossy  problem  having  k  —  k0  —  je,  e  >  0.  We  will  solve  the 
problem  using  Fourier  sine  integral. 

Let 

jo  f°° 

Us(x,X)  =  \-  u(x,y)sm\ydy  (24) 

V  TT  J0 


14 


Then  using  integration  by  parts,  we  see  that 


2  r°° 


J  uyy(x,y)s'm\ydy  =  y -  |u„(x,y)  sin  Ay 

00  ,oo 

—\u(x,  y)  cos  Ay  -  A2  /  tt(x,  y)  sin  \y  dy 

_ r\  ^0 


Because  uy(x,  oo)  — *  0  and  u{x,  oo)  *  0,  we  have 


J  uyy(x,y)  sin  Xudy  =  \J^-Xu(x,  0)  -  X2Us(x,  A) 


-Ay(x)-A2f/S(*,A) 

7 r 


Multiplying  both  sides  of  (23)  with  ^2/*- sin  Ay,  integrating  over  y  =  0  to  oc  and 
making  use  of  (25),  we  have 


-A g(x)  -  A 2Us{x,  A)  =  2jk~Us(x,  A) 

7T 


We  may  rewrite  the  above  equation  as 

~  [[/,(*,  A)c<^»]  =  ±^9(x)e^»  (26) 

Replacing  the  dummy  variable  x  in  (26)  with  r  and  integrating  both  sides  over  r  =  0 
to  x,  we  arrive  at 


t/,(l,A)  =  £7.(0,  A)^'**  +  /^s(t)«<‘W*><~></t 


C/s(0,A) 


[2  f°° 

=  \-  u(0,y)sin  Xydy 

V  7T  2o 

=  \/-  /  /(y)  sin  Ay  dy  =  Fs{ A) 

V  7T  2o 
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■■■  =  (27) 

Finally,  taking  the  inverse  sine  transform  on  (27)  we  get 

[2  roo 

ulx.y)  =  \  ~  /  sin  XyUs(x,X)dX 

V  7T  JO 

=  f°°  Fs(X)e-^2/2jk)xsm\ydX 

V  7T  JA=0 

+  -  —  /°°  A  sin  Ay  f  g^e^^^dr  dX 
^ ir  2jk  Jx-o  yJ r=oyy 

-  -  f  [  /(r)  sin  At  dre~^  ^2:,k^x  sin  \y  dX 

7T  J A=0  J  t=0 

—  [x  g(T)  r  Xsm\ye^/2jk){T-x)dXdT 
Tt2jkJr=0  J  A=0 

=  i  /“  /(t)  r  [cos(r  -  y)\  -  cos(r  +  „)A]  t^'^dX dr 
7 r  Jr—Q  j  A=0 

+i  f  »(T)i  (-|-)  /”  cosXye^^-’UXdr 
7 r  J t—0  jk  \  dy  j  J a=o 

=  if  /(r)  r  [cos(!/  -  r)A  -  cos(y  +  r)A]  e'^'^dX  dr 

7T  Jt=0  •/ A=0 

— 1-  [*  g{r)w~  (  cos  t/A  e-^2|/2'?fc^x-T^c?A|  dr 
njk  J t—0  dy  \  Jo  > 


Defining 


K(x,y,  x0.  y0)  =  /°°  cos(j/  -  yo)Ae~(AW)(x  ro)c?A,  for  x  >  x0, 
Jo 


we  write  the  expression  for  it(x,y)  as 


(x,y)  =  -  I"  f(T)[K(x,y-0,T)-K(x,y]0,-T)}dT 

7T  J  T—0 


4r  /  g{T)-^I<{x,y,T,0)di 

■Kjk  Jt= 0  uy 


njkJr=o*v 'dy  ' 
We  now  evaluate  the  integral  for  Ah  Consider 


/(a)  =  I"  cosaXe-^/2jk)^x-xo)dX,  x> 

J A=0 
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=  r  cosaA 
J  X-0 

=  r  cosaAe"^™ )ftW)l'-i*)d\ 

J\-0 


Because  of  the  exponential  decay,  we  may  differentiate  under  the  integral  sign  to 
obtain 

±1<C)  =  -  r 

da  J a=o 

[e-(A2(x-x0)/2|fc|2)(£-ifc0)] 


(■OO 

J\= 0 


sin 


aA 


d_ 

dx 


(x  —  x0)(e  —  jk0)/\k\2 
sinaAe-A2(r-£o)(e“jfco)/2|fc|2 


d\ 


a\k\" 


A=0 


(x  -  x0)(e  -  jk0) 


(x  -  x0)(e  -  jk0)/\k\2 
[°°  cos  aAe-(A2(l-Io)(£-J,:o))/2i/:!2dA 

Jo 

ja\k\2  f°°  —  X2(x—xn'\/2ik  J\  Jak 


(x  -  x0)k’ 


r  cos  aye-x^x~xo)/2jkd\  =  - 
Jo 

£M  +  ^/(o)  =  o 

da  (x  —  x0) 


(z  -  x0) 


1(a) 


or 


=  0  =t>  1(a)  =  I(a  =  0)e' 


—  \l(a')ejaik^2('x~xo'>'1 
da  1 

Now  /(a  =  0)  can  be  written  as 

I(a  =  0)  =  r  e-{x2{x~Xo)/2W2){i'jko)dX 

Jo 

Let  us  view  this  integral  in  the  complex  A-plane. 

I-Wtk  /W,(A)=  A 


ja2k/(  2(x-xo)) 


-  e 


Rjt^LlprN.  Cff- 

Gmy^A Av^(X)  r  Yifc 

>  R*(A) 


Complex  A-Plane 


(29) 


17 


For  the  integral  to  converge  for  (x  -  x0)  >  0,  i?e[A2(e  -  jk0)}  >  0,  i.e. 

Re[X2(-jk‘)]  >  0  =>  Im(X2k*)  >  0 


or 


or 


or 


0  <  Arg  (A 2k*)  <  7 r 
0  <  2  Arg  (A)  -  Arg  (k)  <  n 


^  Arg  ( k )  <  Arg  (A)  <  -  +  -  Arg  (k) 


Now 


Arg  (k)  =  -  tan  1 

2k0 


kn  ’ 


F<1 


6  A  /\N  *  6 

—  <  Arg(A)  <  -  -  — 


/(a  =  0)  =  /  e->a(*-*o)(-i**)/2|*l2</A 

where  C  lies  in  the  region  of  convergence. 

In  particular,  let  us  choose  a  line  from  0  to  oo  along  the  line  Cw  defined  by 


A  - 


I  2jk 

- 1 

x  —  x0 


w  =  0  to  oo 


(30) 


(31) 


£  _ 
zko 


On  this  path, 


J(a  =  0) 


Arg(A)  =  j+tArg(*0= 


2  jk 


(x  —  Xq)  Jw= 0 


f 

Jw- 


dw 


7 rjk 


2{x  -  sq) 
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Xr,  \  I  r-jk[y-yo)2/[2{x-xo))  132) 

"  x<"yo)  -  f2(x -17) 

Substituting  this  into  (28),  the  field  u{x,y)  for  the  initial-boundary  value  problem  is 
given  by 

«(,.,)  =  i£/ 


1  fx  .  ,  a  /  ti-jA 
-^kLo9iT)dy\i2(x- 


c-iky2/(2(x~T)) \  dr 


It  is  easy  to  see  that 


IT-  =  r  ~TlCOS^-y°)X }e-xH‘-'),mk)d\ 
dx  Jo  2  jk 

54  =  r  -^cosl(y-yoW-^-rW2ik)d\ 

OV2  Jo 


and,  so  for  x  ^  x0 


d2I<  n..8K  n  dl<  _  1  &K_ 

dy 2  1  dx  0r  dx  2 jk  dy2 


which  is  the  same  equation  satisfied  by  u.  Now 


u(x,y) 


27 TX  Jt= 0 


f(T )  6 


.-jk{y-T)2/(2x)  _  ~-jk{y+T)2 /(2x) 


-7TkLs(-T)hK(x'y'  r’0)iT 

=  \l¥~  /”/(T>— 

+  V  27TX  7o  x 

= 

- tt  f  g(T)2jk^-K{x,y,  r,0)  dr 

TTjkJo  OX  0+ 
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in  view  of  (34).  Furthermore,  we  note  that 

d  d 

—K(x,y,xQ,yQ)  =  -~I<(x0y,  x0y0) 
Using  this  and  evaluating  the  first  integral  by  parts,  we  get 


-u(x,y) 


OO  _ 

0  roc 

f(T)e~jkr  /(2x)  -  /  fr(r)t 

Jo 


+  -  I  g(r)-^-K(x,y,  r, 0)  dr 

7T  Jt=0  OT  ,  „+ 


-jkT2/(2x)^7 


V  7TX  l  JO 


-]kT2/(2x)j7 


r  x 

+-  y(r)A'(x,y;T,0)  r_0  -  /  gT(r)K(x,  0;  r,0)dr 

7I'  [  K-0+  ‘/t=° 


=  +v^/(°)  + 


7TX  JO 


fr(T)e-^^2x)dT 


2  /n\  fx  9t ^ 

7T5  V  2a:  V  7T  Jr= 0  \fx  ~T 


=  J—  [/(0) -5(0)] 

V  7TX 

,  /£*[  f°°fT{T):kT2 


/2 jJi  |  r 

7 r  \  Jo 


jkr2 /( 2x 


rz 

^ dr  —  f 

Jo 


$t(t) 


From  the  compatibility  conditions  on  the  initial  and  boundary  values,  we  have 


limu(x,0)  =  y(0)  =  limu(0,y)  =  /(0) 

x— >0  V— >u 


Therefore 


-(x,y)  =  \[^  [4=  I™  fr(r)e-jkT2/{2x)dr-  f  Ldr 

d  n,  V  7 r  Jx  Jt=0  Jr=0y/x-T 

v — ►O'  v 


It  is  easy  to  note  from  (35)  that  for  k  =  k0  —  jt, 


liniic(x>2/)  =0' 
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No  matter  how  small  the  e  is,  we  will  use  the  same  equality  for  the  limiting  case 
of  t  — >  0.  We  will  evaluate  the  integrals  above  approximately  by  replacing  the 
derivatives  with  differences. 


We  now  consider  the  evaluation  of  du/ dy\y—o  at  x  —  2;p-i/2-  Consider  initial  data 
on  the  line  x  =  0.  Let  us  assume  that  this  initial  data  is  known  on  a  uniform  grid 
ym  =  mAy,  m  =  0,1,  •  •  ••  We  approximate  the  derivative  in  the  interval  (ym-i,ym) 
by  the  forward  difference  formula 


df(y) 

~W  ~  Ay  ; 


y  £  (ym-i  >  ym ) 


where  um 
Then 


t(0.ym). 


4=  P  Ur)e-^dr  *  £ 

VX  JT=°  m= 1  A y 

1 

y/x 


1  r  e-jkr2/{2x)dT 
Jr=ym- 1 


Let 
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We  then  have 


y/k/(nx)ym 

_L  [Vm  e-^)/2xdT  =  /  e"Jr"2/2^ 

\/x  J T—ym—l  '  k  * 


y/k/(irx)ym- 1 

y/k/(vx)ym 

\A/(  vx)y  m  —  1 


where  F(//)  =  complex  Fresnel  integral  [8] 


-  OO 

F  /  Mr)e-^^dr  =  £ 

*/T— 0  771=1 


^771  ^771  —  1  \ 

,  Ay 


^  (  \/—  2/m)  -  ^  f  J  —  Vm- 1 
7TX  /  l  V  7TX 


(36) 


Consider  now  the  boundary  data  on  the  line  y  =  0.  Assume  that  we  have  a  non- 

uniform  grid  0,x1,x2)  ’  ‘  ‘  ,xp-i/2  ~  x-  On  the  interval  (xm_j,xm),  we  approximate 

the  derivative  as  . 

.  .  um  -  u771"1 

9t(t)  =  -  t;~t 

771  x771  —  1 

where  u771  =  u(xm,0).  At  the  origin  we  have  u°  =  u0.  We  write 


j: 


-i/2  gT(r ) 


zdr 


=  /""  9^T>  dr  +  P"'' 

/r=0  \/  Xp_i/2  —  r  Jr=xp- j 


-i/2  yT(r) 


=0  yJXp_i/2  X  v^r— 0  -y/Xp-l/2 

We  evaluate  the  first  integral  as 


\JX-p— 1/2  ^ 


zdr 


l 


XP-1  ^r(x) 


^  u771  -  _ 1_ 

TO  — 1  Xyjj  —  i  iim-|  y/%p—\/2 


dr  = 


x  —  r 


:dr 


p— 1  _  x,™-1 


rmr~  (-sv^l37) 

71  =  1  ^  771  —  1 


£m  —  1 


p~r  u771  -  u771"1 


2E 


771  =  1 


XfTl  Xjji  —  i 


^p— 1  /2  X jxi  —  1  Xp— 1/2  ^ 


p-1 


=  2E 


um  —  um  1 


1  y/xp— 1/2  Xm_i  -h  \JXrp— 1/2 


(37) 
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Now 


rxP- 
J  Xp— 1 


1/2  gT(r) 


zdr 


uP-1/2  —  uv~l  [xp-U2  dr 


p-i/^  _  uT~l  rxp_ 

yJXy-\j2  T  Xp- 1/2  Xp- 1  Jx p_i 

uP-l/2  _  UP-1  r; 
{Xp-1/2  *^p— l)  ^0 


%t)— 1/2  *^d— 1  Jxv— 1  \/Xp — 1/2 

yP-i/2  _  [xp-i/2~xp-i  d\ 


Vx 


(p-l/2  _  „P-1 


=  2- 


\JXp—\j2  ~  Xv- 1 

Substituting  (36),  (37),  and  the  above  in  (35),  we  get 

Fp-l/2 


5?/ 


+ 


y =0 


8jfc 


V  7r(Xp_i/2  -  Xp_i) 


u 


xp—l/2 


y=0 


V2 7 


A2/ 


)  )  (uTO  um_i ) 


^  7TXp— 1/2 


2/m 


-  F 


\  'XXp—yj 2 


2/m  — 1 


U  —  li 


m-1 


+ 


*  sj{xv-M2-Xm)  +  ^{xv-\l2-Xm-X)  ^(^"1/2  XP-l) 


8j  A- 


t/'1  (38) 


where 


^p— 1/2 


— (xp_j  +  Xp)  =>  Xp_  1/2  —  Xp_i  —  2^a'P  XP~l)  2  AxP-! 

The  above  equation  is  the  discrete  version  of  a  continuous  boundary  condition  of  the 
form  rs 

+  r(x)u  =  s(x)  (39) 

oy 


where 


s(x)  =  \j2j 


(40) 


8  jk  ^ 


um  _  um-\ 


7T 


and 


^  y/x  —  Xm  +  yjx  —  Xm—\  ^  7r(x  Xp- i) 

F(x)  =  fXe-^2dr 
Jo 


u--1  ,  (41) 


(42) 
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is  the  complex  Fresnel  integral  [8]. 

For  efficient  implementation,  the  PDE  and  the  various  boundary  conditions  will  be 
transformed  into  a  curvilinear  coordinate  system  generated  by  setting  the  lower  ir¬ 
regular  boundary  as  a  constant  curve  curvilinear  coordinate. 
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7.  Transformations  to  a  Curvilinear  Coordinate  System 


Consider  the  narrow  angle  parabolic  equation  with  range  dependent  refractive  index 
(a i  /  0)  given  in  (18) 


— 


— - - T{a;+a2|-  +  |l}»  PDE 

(2jfc0-ai)  l  dy 2J 


together  with  the  tropospheric  boundary  condition 

d 


dy 


u 


(■Em )  2/o)  “t"  r(im)it(xjn )  2/o)  ( ^tti 7  2/o ) 


and  the  impedance  boundary  condition  on  the  irregular  boundary 

+  cju  =  0 


V - 7 - 7 - 7 - 7 - 7 - 7- 


X'YW 

We  will  transform  these  to  a  curvilinear  coordinate  system  (£ ,T] ) 


$  (ju<vvc6*Ay 


C 


_<! _ L 


U. - L 


N  BC  ^ 


t  , - Lr 


C _ d 

— ‘-n 

_ £ _ 

_ / 

i  1 

•7—7— 

- - 7 — 

7 

L  7 — 

yjzo 


Phj'sical  Domain 


Computational  Domain 
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We  assume  that  we  have  a  transformation  of  the  following  form 

z  =  z(f)>  y  =  y(^v) 

The  various  metrics  needed  in  the  transformed  equations  are  [5] 

Shi  =4  +  y\  ,  012  =  x(xv  +  y(yv  . 

We  then  have  with  y/g  —  x^yr,  —  xvy^  =  x^yv  ^  0 

ux  =  -  y^u„)  =  — -  uvy() 

s/9  xiVv 

Uy  =  -^-(-X^  +  X^Ur,)  =  ~ 

v0  yv 


Uyy 


(  uv\ 

i 

Vnn 

— 

—  o 

Urjrj 

\yj 

y2 

y r) 

Substituting  into  the  PDE  we  get 


1 


or 

u(  = 
Letting 


xtVr, 

X(_a\ 


{nv n  ~  unVt)  = 


i _ /  *  .  uv  y^v  Ur)Ti 

(2jko-a1){aiU  2  yn  yl  ”2 


y2 


(  a2  Vim 
-U+ - r- 


xi 


+  —  uv  + 


x( 


(2jk0  -a1)U'T  [V,y„  V*  )  (2Jk o  ~  d)  ’  Vv\  '  (2j&o  -  ai)»? 


&i  = 


xta\ 


(2jk0  -  aj) 


b7  =  -*■ 


f£ 

Vr, 


a  2 


Vnv 


+ 


y2/(2j^o-oi)  xj 


(2jT0  -  fli)y. 


we  express  the  PDE  as 


=  6jU  -f  62^  T  ^3U7j7) 


(43) 


The  normal  derivative  on  a  rj  =  constant  line  is 

012 


uv{  f]  —  const.)  = 


uv - 7===u( 


y/99n 


\!xl  +  y\ 


— - - 11^ - ,  - - 

—  Vtxr)  yx|  +  yjix^yrj  —  y^xv) 


ViVv 


uz 


(44) 


26 


Variation  of  an  arbitrary  vector  r  along  the  rj  =  const,  coordinate  is 

r(_  =  x^x  +  yty  =  a( 


The  unit  vector  along  the  tangent  on  rj  =  const,  is  then 

.  _  jk_  _  +  y& 

y/x$  +  y\  y]x\  +  y\ 


Defining 

xi 

yjx\  +  y\ 

Vi 

\lx\  +  y\ 

(xn  cos  9  +  yv  sin  9)u $ 

Vt,  cos  0  -  Xn  sin  8  cos  9  -  xv  sin  8) 

With  these  substitutions,  the  boundary  condition  at  the  bottom  boundary 
0  gets  transformed  to 

(x„  cos  8  +  y,,  sin  9)  .  m  „ 

Ur,-  — - ,  - +  (j/„  cos  0  -  :r„  sin  0)c!U  =  0  on  rj  = 

V4  + V? 

For  2:^  =  0  and  using  yjx\  +  y\  =  xj  cos  8 ,  we  get 

Un  —  ~  sin  6  cos  8u^  T  Ciy^  cos  8u  =  0  @  y  =  0 


Finally  at  the  top  boundary,  we  have 

+  ru  =  s 


cos  0 


sin  6 


we  express  (44)  as 


ur 


or  | _ _ _ _ 

UT)(t,m,  N)  4-  r(^m,  N)yn{(,m)u{t,m ,  N  )  =  N) 


■v  +  CiU  = 

(45) 


(46) 


(47) 
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8.  Generation  of  Curvilinear  Coordinate  System 

Consider  a  piece-wise  linear  ground  profile  and  a  horizontal  upper  boundary 


Physical  Domain 


Tj-N 


X'rv&pJ- 

•—* - < — 

A/ 

- * - a 

B' 

r* - 45 — 

C ' 

D  oJ-Ol 

LbrvfiV, _ y 

T - 1 

A  Yj 

tl 

■  - —4 

— i 

-7 - r-+- 

- 7 - 

Rectangular, 
Uniform  Mesh 

Computational  Domain 


We  generate  the  (x,  y)  coordinates  of  an  interior  point  by  using  a  linear  interpolation. 
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Letting  Ax,  =  (x,+i  -  xj,  Ay,  =  (y,+i  -  y,),  and  A&  =  (6+1  -  6),  we  write 

Axi 


X  =  Xi  + 


y 


1  - 


Afc 

H 

N 


(*  -  6) 


.  Aj/i. 

»,  +  -  6) 


77 

+  jVSo 


(i<(<  (,-+l 

0  <  T)  <  N 


(48) 


At  any  interior  point,  6  <  £  <  (,+u  the  various  metrics  are  evaluated  as 
A  x; 


xi 


Vi 


A&  ’ 


t) 


(49) 


j  _  V_\  At/,- 


2/t, 


1 

jv 


yo  -  y, 


Ay, 

Afc 


(£  -  6) 


yT-,77  —  o 


N J  A(,  ’ 

At  the  boundary  points,  we  use  the  central  difference  formulas  to  arrive  at 

xt+2  —  Xi  Ax,-+i  +  A  Xi 

W  =  ?«•)  = 


6+2  —  6  A6+i  +  A6 

r  t  _  c  \  _  Ayi+i  +  Ay, 

^  V1  AV  A6+i  +  A<6 


(50) 


(51) 


Note  that  the  analytical  expressions  yield  a  discontinuous  value  for  these  derivatives 
at  the  boundary  points.  We  use  the  analytical  expressions  only  to  generate  grid 
points  and  use  the  central  difference  formulas  to  arrive  at  the  derivatives  w.r.t.  (. 
In  other  words,  once  the  grid  points  are  generated  on  the  lines  AA',  BB\  . . . ,  etc., 
we  assume  that  the  space  is  smoothly  connected  through  the  grid  points.  In  the 
numerical  implementation  using  Crank-Nicolson  implicit  scheme  [6],  the  metrics  are 
needed  at  the  midpoint  w.r.t.  £,  i.e.,  at  £  =  6  +  A£,-/2  and  the  interior  point  formulas 
are  applicable.  For  a  uniform  mesh  in  the  computational  domain,  A 6  =  1,  rj  —  q, 
</  =  0. 1, 2,  •  •  • ,  Ar 

Ax,  (52) 

(53) 

(54) 


H  = 
Vi  = 

Vn 

V  = 


1  -  Vi A*'  ^ 

1  (  y,+i  +  y. 
y  o  - 


y«(?  + !)  -  vdq)  = 


Ay, 

’  N 


N 


y»+ 1  +  y; 


2 

y«+i  +  yi 


)  +  7/y° 


+  qyn  =  yo  +  {q-  N)y„ 


(55) 


so  that  y(q  +  1)  -  y{q)  =  yv- 
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9,  Numerical  Implementation  by  Crank-Nicolson  Scheme 

Consider  the  narrow  angle  PE  together  with  the  boundary  condition 

u ^  =  b\u  -f-  i>2 )  T  (PDE) 

un  _  sin  6 cos  0J  +  c^cos  6u  =  0  at  rj  =  0  (BC2) 

uv  +  ry^u  =  y-r,  s  at  rj  =  N  (BCi) 

We  would  like  to  implement  the  above  using  a  Crank-Nicolson  implicit  scheme  [6] 


At  -  1 


We  use  the  notation  uvq  =  u(£p.rjq)  =  u(xp.yq) 

The  various  derivatives  assuming  =  Arj  =  1  are 


U^p-l/2;T)q)  = 


<  ~  upq  1 


—  u-n{ip-\/2-,  Vq)  ~  77  [Ur,(£p-1>  Pq)  +  UV  {(p’Vq)] 


1  r  p— i  „,v~  i  1  „,p 


‘9+1 


-  <_i  +  <+i  -  Vi 


or 


u 


(p~  b9) =  1 1“’+' +“'+i  _  +u«_i) 

1  r<+1  -  2<  +  +  <;!  -  2„;-'  +  <:! 


m  \P  2'^)  2 
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Substituting  into  the  PDE,  we  get 


<  -  <' 


,  ,  i  \ + “r1 
2 - 


+ 


62  {p~\'q;  t 

7 - -  (<+1  +  u 


p  _  ,,p 

9+l  U9 


63  p  (uj+1  _  2u;  +  +  <;{  -  2 u\-'  +  <:]) 


Rearranging  the  terms  we  get 


\  (y  +  ^3)  U5+1  _  f1  +  63  ~~  2")  uPq  +  2  _  ^  W®_1 


bo 


+  ^  (y  +  &3^  <+l  +  f1  -  ^3  +  7-)  <  1  +  O  f^3  -  -7)  <-l  -  0 


Now  let 


V 

2 


1  ft  1  ,\p_1/2 

a  =  2O+2M, 


*  =  (+-| 


P-1/2 


IL  -b{ 

2  l  3  2 


P-1/2 


We  then  have  for  p  =  1 , 2,  •  •  •,  q  =  1 , 2,  •  •  •,  A  —  1 

«<+i  -  (i  +  +  7<-i  =  -  (!  -  /?K_1  -  (56) 

However,  we  will  extend  the  applicability  of  this  equation  over  q  =  0,1,-  --N  to 
accomodate  the  derivative  boundary  conditions  on  the  lower  and  upper  boundaries. 
For  q  —  0,  we  have 


aup  -  (1  +  (3)up0  +  7u!_!  =  -aup  1  -  (1  -  P)upQ  1  -  7^ 
for  q  =  N,  we  have 


p-i 

1 


aupN+1 


(1  +  P)uPN  +  7uN-1  =  ~auN+ 1  _  (!  ~  P)UN  1  -  luN-l 


,p-l 


(57) 


(58) 
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From  the  BC2  at  the  lower  boundary  we  have 


Multiplying  this  with  4(7)0  ^  an<^  adding  to  (57) 


(a  +  7)1*1  -  ^1  +  P  -  27  Ci yri  cos  0  2^-  sin  6  cos  6 


From  the  BCi  at  the  top  boundary  we  have 


(59) 
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1  f  ^+1  ~  ^JV-1  mN+1  u7V-1 


+ 


16;  k  lix  +  u1^1 


M  7tAxp. 


J/rj 


I  <r^  /  \ 

—  2/p  {  ~X  L(Um  Um-l) 

m= 1 


^  7TXp_i/2 


Urn  F 


^  7 rXp—i/2 


Vm  —  1 


m  _  ,.m-l 
“N 


^  m=l  (\/2:p-l/2  —  +  y/xp-\j2  ~~  xm-lj  \|  7r^XP~ 


,  ,  1Qik  P-1 

+  A  I  ~ 7  ‘  UN 


or 


w+\  +  <+i  -  (<-i  +  <-i)  +  ^  (u*  +  u*  *)  -  yp4,sP  1/2 


multiplying  with  —  1|/2  and  adding  to  (58) 

1  +  P  +  yv8a 


jk 


UN  +  (7  +  aOuA’-i 


1  -  P  ~  Vv 8a, 


jk 


\  ttA^p-j 

~  (7  +  ~  4cn/77sp-1/2 


Letting 


^  7T AXp 


A  =  /9  +  ?/p8a, 


jk 


^  7T2^iTp  —  5 


7'  =  7  +  a 

we  rewrite  the  above  equation  as 


(1  +  A )upN  -  7'ua--i  =  (1  -  ^)un  1  +  T^jv-i  +  4q2M 


P-1/2 


where 

0P-1/2  - 


V2j 


Ay 


)  ")  (^m  Hm  — l) 


/ 

k  \  ( 

^s. 

3 

1 

-» 

v\ 

TTZp-l/2  /  \\ 

^  ’KXp—\j2 


Vm  —  1 


'8;^  ^ 


UAr 


^  m= 1  xp-\j2  -I'm  "h  y/xp- 1/2  xm-l'j 


16  jk 


^  7rAx 


p-i 


p-i 


(60) 


33 


and 


Xxp— i  —  %p  %p—  1 

Xp_  1/2  =  Xv-\  =  ~{xp  +  Xv-\) 


We  augment  the  equations  given  in  (56)  for  g  =  1,2, ...  ,7V-  1  with  (59)  and  (60)  for 
q  =  0  and  Ar,  respectively  to  define  for  9  =  0, 1, 2, ...  ,7V.  The  system  of  equations  so 
defined  can  be  expressed  as  a  matrix  equation  of  the  form 


'XX 

r  v 

Uo 

V 

«1 

XXX 

XXX 

* 

><: 

_ 1 

>■«  ••• 

‘XX 

r  p—1  1 

u0 

XXX 

p  —  1 
< 

XXX 

X  X  _ 

- 1 

••• 

S3 

1 

(61) 


4  ayvsp  1/2 


where  X  denotes  a  non-zero  entry.  The  tridiagonal  matrix  on  the  left  hand  side  of 
(61)  can  be  inverted  efficiently  to  yield  a  solution  on  line  £  =  in  terms  of  the  field 
values  on  the  line  £  =  £p-i-  Equation  (61)  can  be  used  to  march  forward  in  range 
starting  from  initial  data  specified  on  £  =  0. 
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